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We investigate the resilience of symmetry-protected topological edge states at the boundaries 
of Kitaev chains in the presence of a bath which explicitly introduces symmetry-breaking terms. 

Specifically, we focus on single-particle losses and gains, violating the protecting parity symmetry, 
which could generically occur in realistic scenarios. For homogeneous systems, we show that the 
Majorana mode decays exponentially fast. By the inclusion of strong disorder, where the closed 
system enters a many-body localized phase, we find that the Majorana mode can be stabilized 
substantially. The decay of the Majorana converts into a stretched exponential form for particle 
losses or gains occuring in the bulk. In particular, for pure loss dynamics we find a universal exponent 
a ~ 2/3. We show that this holds both in the Anderson and many-body localized regimes. Our 
results thus provide a first step to stabilize edge states even in the presence of symmetry-breaking 
environments. 

PACS numbers: 03.65.Yz, 67.85.Lm, 71.10.Pm 

I. INTRODUCTION 

Symmetry-protected topological states are a well es¬ 
tablished class of phases of matter which encompasses 
a rich variety of microscopic incarnations, ranging from 
spin chains to fermionic systemJIH^. These phases of mat¬ 
ter are a direct manifestation of a (non-local) symmetry 
in a model Hamiltonian, first identified in the context 
of the Haldane phase of spin-1 chainwhich can give 
rise to topological order and zero energy modes at the 
sample boundaries. A prominent example of this mech¬ 
anism is provided by the Kitaev chairP which exhibits 
Majorana edge stateP once defined on a open geometry. 

These edge modes, which have anyonic statistics and can 
be potentially useful for quantum information purposed, 
have a direct relation with a global Z 2 parity symme¬ 
try of certain one-dimensional superconductorP, usually 
granted by the proximity effect. However, naturally oc¬ 
curring dissipation can explicitly break the protecting mi¬ 
croscopic symmetr^EEl. In this context the major chal¬ 
lenge is to stabilize the Majorana modes by exploiting 
additional mechanisms which could counteract the dras¬ 
tic action of symmetry breaking terms in open quantum 
system settings. 

In this work, we investigate the resilience of Majo¬ 
rana edge modes in the presence of incoherent symmetry¬ 
breaking dissipation. We show that strong disorder can 
lead to an exponential gain in the edge mode’s stability 
compared to the case of a homogeneous system. While 
in a homogeneous system we show that the temporal de¬ 
cay is exponential, in presence of disorder we find that 
the Majorana mode occupation decays in a stretched ex¬ 
ponential form. Remarkably, when the coupling to the 
environment induces pure particle losses, the associated 
exponent a 2/3 is universal and independent of any 
microscopic details within our numerical accuracy. 


The effect of disorder on the time-dependent dynam¬ 
ics of closed quantum many-body systems has been ac¬ 
tively investigated in both non-interacting and interact¬ 
ing scenarios. For noninteracting particles in low di¬ 
mensions, disorder leads to Anderson localization (AL): 
all single-particle eigenstates becom e localized immedi¬ 
ately for infinitesimal disordei^^^^^. Recently, the con¬ 
cept of AL has been generalized to interacting systems in 
the context of many-body localization (MBL)P“i^. Im¬ 
portantly, compared to quantum many-body systems in 
thermodynamic phases, MBL systems display unconven¬ 
tional propertiesSHm^ particular, this encompasses 
the possibility of phas e tran sitions at nonzero temper¬ 
ature in one dimensiorP^^I^H^. Moreover, it has been 
shown that topological properties, in one-dimensional 
thermodynamic systems only present in ground states, 
can survive over all of the spectrunP^^^^, greatly enhanc¬ 
ing the reg ime of stability for topological order in closed 
systems^. 

Motivated by these findings, we address in this work 
the fundamental question if and how the MBL mecha¬ 
nism can affect the stability of topological order, and in 
particular the resilience of edge modes, when the sys¬ 
tem is coupled to an incoherent bath. This question 
finds natural applications in diverse systems, such as solid 
state devices and cold atom settings, which have inherent 
sources of dissipation. Our results support the conclusion 
that MBL qualitatively enhances the resilience of edge 
states in symmetry-protected topological states even in 
the presence of the worst-case scenario of a bath inducing 
symmetry-breaking perturbations. 

In concrete, we study the combined effects of disorder 
and interactions in a Kitaev chain coupled to a Marko¬ 
vian bath which explicitly breaks the parity symmetry of 
the system. The setup we have in mind is portrayed in 
Fig. Ill and is described in detail in Sec. |n| Our analysis 
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focuses on the decay of the Majorana edge mode under 
the effect of particle losses, starting from an initial state 
where such mode is occupied. In Sec. Ill, we first of 
all study the case of a homogeneous system as a bench¬ 
mark, where we show that the Majorana mode decays 
exponentially as a function of time. Specifically, we con¬ 
sider a particular setup where the dissipation acts only 
on a part of the system (red area in Fig. with a vary¬ 
ing region near the boundaries that is not coupled to the 
Markovian bath. This represents an idealized model sys¬ 
tem which, however, is potentially relevant for a variety 
of experimental systems as we discuss in detail in Sec.[n| 
Both in the presence and in the absence of interactions, 
we find that for the case of particle losses the survival 
time of the Majorana edge mode increases exponentially 
with the number of protected sites, that is, the number 
of sites close to the edge where no dissipation occurs. 
This increased lifetime can be immediately traced back 
to the nature of the edge mode wave function which only 
exhibits an exponential tail extending into the bulk. 

In Sec. |IV| and Sec. |V| we discuss how disorder dras¬ 
tically affects the above mentioned scenario. As antic¬ 
ipated before, in contrast to the exponential decay in 
the homogeneous system, the Majorana edge mode de¬ 
cays as a stretched exponential with a universal exponent 
a 2/3 for the pure loss dynamics. This points to¬ 
wards an increased robustness of the edge mode when 
compared to the clean case. In particular, this universal 
behavior occurs both in the Anderson and many-body 
localized regimes, i.e., also in the presence of interactions 
which makes the system generic. This shows that dis¬ 
order, and, in general, localization effects, which induce 
non-ergodic dynamics in the system, are not detrimental 
to the survival of Majorana modes in the presence of a 
bath, and on the contrary, greatly enhance their stability. 
This identifies disorder as a key mechanism to stabilize 
topological features in open quantum systems, and thus 
to systematically improve the stability of potential quan¬ 
tum memories even in the presence of dissipation which 
explicitly breaks the symmetry protecting topological or¬ 
der. Finally, in Sec. |V| we present a case study for the 
loss and gain case, where a qualitatively similar behavior 
can be observed as well. We find that the decay of the 
Majorana again is of stretched exponential form. 


II. MODEL HAMILTONIAN AND MASTER 
EQUATION 

The Kitaev chain can be realized in a variety of differ¬ 
ent experimental architectu res including AMO as well as 
solid-state setups Hamiltonian describes 

spinless fermions on a lattice of N sites with a pairing 
ternP: 
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Figure 1. (color online) Schematic picture of the dynamics 
of a disordered Kitaev chain subject to a bath. The system 
is initialized in the ground state of the ideal Kitaev chain, 
where it is better visualized splitting the fermionic degree of 
freedom at each site into a pair of Majorana fermions (small 
full circles). In the ground state, all Majorana fermions in the 
middle couple over bonds, leaving the first and last Majorana 
fermion uncoupled (blue circled). The dynamics is the subject 
to a disordered potential (black, thick line) in addition to 
pairing to a superconductor (grey shaded area) and to the 
action of a bath acting on the bulk (red area). The presence 
of the bath explicitly breaks the parity symmetry of the wire 
due to particle loss. 


where cj,ci are fermionic creation/annihilation opera¬ 
tors at site I = The first two terms repre¬ 

sent tunneling between lattice sites as well as pairing 
induced by a superconductor, whose magnitude can, in 
principle, depend on the specific lattice site, and the last 
term is a chemical potential. In the homogeneous regime 
Ji = J, Ai = A, |/i/J| < 1,A 7 ^ 0, the model hosts 
a symmetry-protected topological phase directly related 
to the exact parity symmetry of the Hamiltonian, that 
is, [H,P] = 0,p = nili(l “ 2n;), where n; = cjci is 
the local particle number operator. In order to assure 
that the properties of the system are generic we account 
for a weak, parity-preserving nearest-neighbor interac¬ 
tion making the system non-integrabl^^. This yields as 
the full Hamiltonian: 
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with U/J 1. Notice that this model can be mapped 
via a Jordan-Wigner transformation to a transverse-field 
XYZ chain, and as such is nonintegrable and genericP^. 
In the case where ja = U = 0, J = A, the ground state 
(GS) of the open Kitaev chain supports a localized Ma¬ 
jorana edge mode (see illustration in Fig. whose wave 
function is restricted to the edge lattice sites onl}!^. In 
the odd parity sector for odd number of sites, the GS 
reads 


IV^o) = 2-^ (|1), + |0)J + (|1), - |0),)] . 

( 3 ) 

where the occupation of the zero energy mode, the edge- 
edge correlation, 0 satisfies: 

6» = (V»o|0|^/’o) = 1 , 0 = i (ci + 4) (cat - ct) ■ (4) 
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and 0 denotes the Majorana mode occupation operator. 
If we would slightly perturb our Hamiltonian with > 
0 , but still sufficiently small not to leave the topological 
phase, the Majorana zero mode develops an exponential 
tail extending into the bulk, with mo st of i ts weight still 
located at the boundaries of the chaiiP^^^. 

The above particular choice of |'0o) represents an ideal¬ 
ized initial condition which, of course, relies on an exact 
implementation of the Hamiltonian, or on a perfect dis¬ 
sipative state preparation of the state of interest (along 
the lines of Ref. [36]) . As we would like to argue at the end 
of this section, however, our full setup is still sufficiently 
general in order to describe the generic dynamics. 

Remarkably, it has been recently shown how the topo¬ 
logical features of the ground state in one-dimensional 
systems - in particular, the presence of edge modes and 
a degenerate entanglement spectrum - can be extended 
to the entire spectrum in the presence of strong disor- 
(^e] i2i | 23 | 2 4i^^ This mechanism, a direct consequence of 
many-body localization, can lead to a stabilization of 
topological features at nonzero energy density above the 
ground state in a closed system. 

Our main focus here is to investigate the fate of the 
Majorana zero mode in the presence of symmetry- break¬ 
ing dissipation. The microscopic details of the particle 
loss and/or gain mechanism will strongly depend on the 
experimental architecture. In a cold atom experiment, 
losses are naturally occurring as either inelastic scatter¬ 
ing processes due to the inime rsion in a BEG of molecules, 
or background collisionJ^lEII. These effects are often very 
well modeled as a Markovian batlP^. In solid-state sys¬ 
tems, loss and gain terms can emerge as a consequence of 
a coupling to a grand-canonical bath. For Kitaev chains 
realized by placing carbon nanotube wires on top of a 
p-wave superconductor, particle losses and gains can oc¬ 
cur by electron tunneling into and from the substrate. In 
this work, we concentrate on the scenario where both the 
loss and gain can be described within a Master equation 
framework. At a qualitative level, this implies that the 
parity symmetry in the system is violated at all energy 
seales - which we interpret as a worst-case scenario for 
symmetry-protected topological states. The system dy¬ 
namics can then be described using a Lindblad master 
equation of the form: 
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with p the density matrix of the system, ni and 7 / are the 
local loss and gain rates, respectively. We consider uni¬ 
form dissipation, but restricted to the bulk of the chain 
where = /i:, 7/ = 7 for d < / < A" — d, and ki = 'yi = 0 
otherwise with d the distance of the first dissipative site 
from the edge. For an illustration, see Fig. In other 
words, we mainly consider the influence of dissipation in 
the bulk, not at the boundaries. 


The main motivation for the choice of the setup is 
twofold. First, from a theoretical point of view, by re¬ 
stricting dissipation to the bulk we address the funda¬ 
mental question of whether information stored in Ma¬ 
jorana modes is stable against symmetry-breaking pro¬ 
cesses in the bulk. While dissipation acting on the edges 
of the system opens a direct channel for the decay of 
the Majorana mode, protecting the boundary sites al¬ 
lows us to extract the bulk contributions and therefore 
to unravel the influence of the indirect channel. Second, 
our setup might be relevant also in a general experimen¬ 
tal context where dissipation is not uniform and particle 
losses or gains are not occuring homogeneous throughout 
space. In that case, isolated regions with weak or even 
approximately vanishing decay rates can form for which 
our setup represents an idealized reduced model system. 

We remark that we have checked via extensive numer¬ 
ical simulations whether the specific choice of the distri¬ 
bution of the decaying sites can have an impact on the 
dynamics and our main results. We found that solely 
the number of sites which are not affected by dissipation 
mattered in terms of the qualitative behavior - as long as 
the edge sites are not included. 

The figure of merit in our study is the fate of the Ma¬ 
jorana mode: 


0{t) = tr:[p{t)e], (6) 

as a function of time t in presence of the incoherent 
bath. For our initial density matrix p{t = 0) = \ipo){ipo\ 
with ^ = 1 , we numerically solve the master equation 
in Eq. For the noninteracting system with = 0, 
the resulting master equation is of quadratic fermionic 
form a nd ca n therefore be solved efficiently for large 
systemJ^niini^ We use this noninteracting limit as a bench¬ 
mark for our more general and exact numerical imple¬ 
mentation where interactions can also be incorporated 
at the expense of limited system sizes. We find, however, 
that, for the considered scenarios, finite-size effects can 
be neglected (see discussion below). We start by com¬ 
puting the dynamics of 0{t) in the clean case, and then 
discuss the role of disorder, showing how the latter case 
displays to an exponential gain in stability of the Majo¬ 
rana mode when compare to the former one. 

Let us note that 0 in Eq. ^ does not measure the 
overlap with the idealized wave function of the Majorana 
mode with respect to the final Hamiltonian for the gen¬ 
eral parameter set considered in this work (the wave func¬ 
tion will be generically localized at the edge, with a finite 
contribution in the bulk), but rather those contributions 
which are located right at the edges. This, however, is 
not a shortcoming as we will now shortly discuss. First 
of all, the dynamics at the edges has to reflect the dy¬ 
namics of the full Majorana wave function because most 
of its weight is located right there. Moreover, we expect 
that the considered scenario is also the experimentally 
most relevant one. In particular, measuring operators lo¬ 
cated on one lattice site is much easier than measuring 
an extended object such as the full Majorana wave func- 
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Figure 2. Impact of particle loss on the edge-edge correlation 
0 in the clean case = 1: (a) with a finite external field 
H = 0.02, 0.03, 0.04 {U — 0) and = 1/2 for one protected site 
d—1] (b) with an on-site repulsion U = 0.04,0.05,0.06 (/x = 
0) and = 1 for d = 1. Circles and boxes in (a) and (b) show 
the dynamics for systems of different lengths (N=9,ll). 


tion. In this context, it is also important to emphasize 
that the precise form of the Majorana wave function in 
a concrete experiment is not a priori known because it 
depends on a variety of experimental details - especially 
in the presence of interactions. 

Finally, these arguments also motivate our specific 
choice of initial condition in Eq. In any realistic 

experimental system it will, of course, not be possible 
to precisely prepare the considered state due to experi¬ 
mental imperfections. As long as the Majorana mode is 
initially occupied, however, our observable 0 is capable 
to detect it without knowledge about microscopic details. 
Therefore, we expect that these details only lead to quan¬ 
titative and not qualitative changes. 


III. HOMOGENEOUS SYSTEM 


Using the model system described in the previous 
section we now aim at studying first the influence of 
symmetry-breaking particle losses for the homogeneous 
system. This will serve as a benchmark for the further 


analysis in Sec. |IV| where we will show that the Majorana 
edge mode can be stabilized by the inclusion of disorder. 
For the homogeneous system we analyze the noninteract¬ 
ing and interacting cases separately in order to single out 
the potential influence of integrability-breaking interac¬ 
tions. 

The main result for the homogeneous system is that 



Figure 3. Panels (c -d): exponential dependence of the effec¬ 
tive decay time F on the number of protected sites d, with N 
= 43 (11) in the left (right) panel. 


the Majorana mode occupation decays exponentially in 
time in consequence of the breaking of the protecting 
parity symmetry. However, we find that, remarkably, the 
associated decay rate F is highly sensitive to the distance 
d of the first dissipative site from the edge, but not to 
the microscopic parameters. In particular, interactions 
do not lead to a qualitative change of the dynamics but 
rather only increase the decay rates. Our main findings 
are summarized in Fig. 


A. Non-interacting case 

We start by considering the benchmark of noninter¬ 
acting particles with U = 0. In this limit we have used 
two different methods to solve the Master equation. On 
the one hand we have made use of the property that the 
Master equation at U = 0 becomes a quadratic form in 
terms of fermionic operators. This allows us to use meth¬ 
ods developed in Ref. Eg such that we are able to obtain 
the solution for large lattices up to = 100 sites. On 
the other hand we have solved the Master equation di¬ 
rectly on the basis of a Runge-Kutta algorithm for up 
to A" = 11 lattice sites which can also be extended to 
interacting systems studied later on. By comparing the 
results of both methods, we have found that finite-size 
effects are negligible in the parameter regime considered 
here. 

As anticipated before, we find that in the non¬ 
interacting homogeneous system the Majorana mode oc¬ 
cupation 0{t) decays exponentially in time: 

(7) 

with F the decay rate, both in the cases of a global bath 
acting on all sites, ki = /^, and in the case of a bath which 
does not act on a small region close to the edge, ni = n 
for / G [d + 1,1/ — d] (with d the number of protected 
sites close to the edge), and 0 otherwise. In the follow¬ 
ing, we focus on the bulk contribution, investigating the 
Majorana decay as a function of d. 

In Fig. Hf, time traces of the Majorana occupations 
are shown for the case where the first dissipating lattice 
site is located right next to the edge, i.e., d = 1. As one 
can see, the decay rate depends crucially on the strength 
of the chemical potential. This can be traced back to a 
fundamental property of Kitaev chains at /x = 0. There 
the Majorana qubit decouples completely from the bulk 
chain and Majorana mode occupation operator 0 com¬ 
mutes with the Hamiltonian =0. As a conse¬ 

quence, we have that 0{t) = 1 for all times t when /x = 0. 
Therefore, at this particular point in parameter space 
the Majorana completely decouples from the dissipative 
dynamics and becomes completely stable. However, this 
stabilization mechanism depends crucially on parameter 
fine-tuning and any perturbation, always present in a re¬ 
alistic experiment, induces an exponential decay. 

The infinitesimal sensitivity of the Majorana qubit sta¬ 
bility can be seen by considering the effects of a non- 
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vanishing chemical potential /i > 0, which is shown in 
Fig-i Then, the Majorana qubit experiences an effec¬ 
tive coupling to the bulk modes and the exponential de¬ 
cay induced by the symmetry-breaking losses sets in. The 
decay rate grows with increasing chemical potential. 

Remarkably, it is possible to increase the lifetime of 
the Majorana mode substantially also for the homoge¬ 
neous system without parameter fine-tuning. Specifi¬ 
cally, we find that the decay rate is exponentially sen¬ 
sitive to the distance d of the edge from the first dis¬ 
sipative lattice site. This is shown in Fig. where T 
is plotted as a function of d. The decay rate is ex¬ 
tracted numerically from two time points t 2 > ti with 
T = Log[0{t2)/0{ti)]/{ti — ^ 2 ). We have checked that 
the resulting T’s are independent of the precise choice of 
ti and t 2 as long as we are in the asymptotic long-time 
regime. As a consequence, in the case of pure particle 
losses, the bulk only contributes weakly to the decay of 
the Majorana. Notice that this is different from the case 
with losses and gain discussed in Sec. |V| We attribute 
this strong sensitivity of the decay to the exponentially 
small overlap of the Majorana wave function with the 
loss mechanism on distant lattice sites. In particular, as 
long as the evolved Hamiltonian stays in the topological 
phase, the Majorana wave functions are exponentially lo¬ 
calized in the vicinity of the edges. Therefore, the direct 
influence of a particle loss at a given lattice site on the 
Majorana mode is exponentially small. However, excita¬ 
tions that are created by a particle loss could, in prin¬ 
ciple, still propagate to the edges potentially inducing a 
decay of coherence. The numerical data, however, sug¬ 
gests that the dominant principle is solely the overlap 
of the Majorana wave function on the first dissipative 
lattice site. 


B. Interacting case 


After having discussed in detail the noninteracting ho¬ 
mogeneous system we now aim at analyzing the influence 
of interactions. This is important for clarifying which of 
the effects are generic and independent of particular pa¬ 
rameter choices. The numerical results for the interact¬ 
ing case are illustrated in Fig.|^ where, along the lines of 
the previous section, we check the resilience of the edge 
modes for increasing interaction strength. 

As for the noninteracting case, we find that the Ma¬ 
jorana edge modes decay exponentially in time, see the 
numerical data presented in Fig. Again, the decay 
rate decreases exponentially with the distance d of the 
first dissipating lattice site, see Fig. In contrast to 
the noninteracting case, however, an even-odd effect of 
F as a function of d is visible. Importantly, the model 
is now generic which becomes apparent at the fine-tuned 
parameter point /r = 0 of the noninteracting model, see 
Sec. [Thai where the Majorana mode is stable for d > 1. 
In Fig. we show time traces of 0(t) at /i = 0 but non¬ 
vanishing interactions. As one can see, interactions now 


lead to a decay of the edge mode. 

Before we study the impact of disorder, we summarize 
the main result of this section: Symmetry-breaking per¬ 
turbations lead to a decay of the Majorana qubit even 
if the edge states are protected from particle loss, in¬ 
dependently on the presence or absence of interactions. 
However, the survival time of the edge mode increases 
exponentially with the number of protected sites d. 


IV. MANY-BODY LOCALIZED KITAEV CHAIN 

For closed systems, it has been recently observed that 
disorder provides a generic mechanism to stabilize topo¬ 
logical order throughout the f ull many -body spectrum, 
i.e., up to infinite temperatnr d^^ * ^^ ^ Strikingly, this 
holds even for one-dimensional systems where long-range 
(topological) order at nonzero temp erature is not pos¬ 
sible in thermodynamic phases^^^^^^^. These many- 
body localized systems are non-ergodic and consequently 
promise a key mechanism to stabilize qubits in closed 
quantum system settingd^. 

For Anderson-localized systems it is known that a 
coupling to low-temperature baths in general induces 
a nonzero conductance in the context of variable-range 
hoppin^^. However, depending on the spectral details 
of the bath, localization can persist in the sense that the 
coupling between system and bath becomes irrelevant at 
low energied^. For small baths, not of thermodynamic 
nature, it has been show n tha t localization can persist 
in a strong disorder hmi4^^^. In many-body localized 
systems the influence of low-temperature thermal baths 
onto the broadening of local spectra has been studied re- 
centl}^^. The static properties of MBL systems coupled 
to small baths and the transition of the combined sys¬ 
tem to a Wigner-Dyson statistics have been addressed in 
Ref. im 

In the following, we study the many-body localized 
Kitaev model in presence of an incoherent structure¬ 
less bath described by the Master equation in Eq. 

We show the stabilization property of many-body local¬ 
ized phases remarkably extends also to this worst-case 
open quantum system scenario, thus identifying non- 
ergodicity of closed systems as a possible route to stabi¬ 
lize edge states in the presence of detrimental (symmetry¬ 
breaking) noise. Our main result, summarized in Fig. 
is that the Majorana mode can be stabilized by the in¬ 
clusion of disorder with an exponential gain in resilience. 
Specifically, the exponential decay of the homogeneous 
system converts into a stretched exponential decay in 
presence of disorder, that is: 

{0{t)Us ( 8 ) 

with (... )dis denoting the average over disorder realiza¬ 
tions. We find that within our numerical accuracy the 
associated exponent a is universal and close to 

a - 2/3. (9) 
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In our numerical implementation, we choose bond dis¬ 
order for simplicity by introducing inhomogeneous cou¬ 
plings Ji G [0,11^], I = 1... N — 1 sampled from inde¬ 
pendent uniform distributions of width W = eJ. We 
will characterize the strength of the disorder via J dif¬ 
fering from W by Euler’s number e such that the critical 
point of the noninte ractin g model is located at ji/J = 1 
as in the clean cas d^^ l ^^ i In the following simulations, 
we are interested in the large disorder limit such that 
the localization length of most eigenstates is limited to 
few lattice spacings (thus limiting finite-size effects). For 
each parameter regime, the results are averaged over up 
to 10 000 disorder realizations in order to ensure conver¬ 
gence. We expect, however, that the generic features are 
independent of this precise choice of bond disorder, in 
particular, because it is known that the properties of the 
bond-disordered Kitaev and the related Ising chain are 
equivalent to the case of site disorder. 

Before diving into a detailed discussion of our results, 
let us comment on the case where the dissipative losses 
act on every lattice site, i<ii = k, in particular, includ¬ 
ing the edges where the Majorana modes are located. In 
this case, we find that the Majorana mode occupation 
still decays exponentially as found for the homogeneous 
case. Thus, the Majorana qubit is not stable against in¬ 
coherent onsite losses as one might expect. This local 
decoherence mechanism cannot be stabilized using disor¬ 
der. However, as we will show below, protecting sites in 
the immediate vicinity of the boundaries lead to a sub¬ 
stantial stabilization. 

We start by discussing as a warm-up the noninteracting 
(Anderson) limit with nonzero chemical potential /i > 0 
but U = 0. Our numerical results on this case are shown 
in Fig.[^. In panels[^-d, we plot the same quantity for a 
clean (c) and disordered system (d, same data as in a) as a 
function of the rescaled time (Jt)^/^. While in the clean 
system the decay stays exponential as discussed in the 
previous section (the curve bends down in lin-log scale), 
in the disordered case one gets an almost flat line (but 
for additional coherent oscillations), demonstrating an 
exponentially increased stability of the edge mode, and 
a decay according to a 2/3. The decay becomes faster 
for larger chemical potential as this increases the coupling 
between the Majorana qubit and the bulk modes. An 
analysis of the time derivatives of the Majorana mode 
survival also confirms this decay. 

The same stabilization mechanism emerges in the 
many-body localized regime, where the Majorana oper¬ 
ator is expected to be renormalized, still be localized at 
the boundar}EI. The resulting decay of the Majorana 
mode for different interaction strength for the case fi = 0 
is illustrated in Fig. S', f, while panels (e-d) show the 
comparison between clean and disordered case. Again, 
the decay is of the form of Eq. ^ with a ~ 2/3. 

Finally, let us consider the case of both nonzero in¬ 
teractions and nonzero chemical potentials. The result¬ 
ing temporal behavior of the Majorana mode is shown 
in Fig. Here we consider a comparable strength for 





Figure 4. Impact of particle loss on the edge-edge correlation 0 
in the disordered case Ji G [0, e]: (a) with a finite external field 
ji = 0.02, 0.03, 0.04 {U = 0) and = 1/2 for one protected site 
d = 1; (b) with an on-site repulsion U — 0.04,0.05,0.06 (/x = 
0) and = 1 for d = 1. Circles and boxes in (a) and (b) 
show the dynamics for systems of different lengths (N=9,ll). 
Panels (c-f) show the time trace with a scaled time-axis to 
illustrate dependence. In panels (c) and (e), we show the 
clean case, where the decay is exponential, and thus, in the 
rescaled axis, an algebraic decay with power larger than 1 is 
observed. In panels (d) and (f), we contrast this with the 
disordered case, where the decay is a stretched exponential 
with a power a ~ 2/3. 


both U = J/5 and /x = J/5 at a dissipation strength 
n/J = 3/2. For comparison also the respective homoge¬ 
neous case is shown that exhibits the pure exponential 
decay already discussed in Sec. |HI[ The many-body lo¬ 
calized case on the other hand again shows the stretched 
exponential decay, with an additional coherent oscillation 
on top of the decay is found. 


V. EFFECTS OF PARTICLE GAIN 

So far, we concentrated on a setup where only parti¬ 
cle losses occur. Such as scenario corresponds mainly 
to potential implementations based on cold atom se¬ 
tups, where the main dissipation channel is given by 
inelastic scattering between molecules in the reservoir, 
and particles in the wire. In a condensed-matter set¬ 
ting, a more natural environment would generically be 
a grand-canonical bath inducing both particle losses and 
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Figure 5. Limit of strong interaction and external field fi = 
U = 0.2 with d=2 protected sites and k = 1.5 for the clean 
and disordered case. Here, disorder increases the oscillation 
amplitude and slows down in order of magnitude the decay of 
the edge-edge correlation. 

gain^^^l^. This is the scenario we focus on in this sec¬ 
tion. Specifically, we consider equal rates ki = 7 / for 
both losses and gains for simplicity. Of course, baths in 
typical condensed-matter systems will not be exactly of 
the type described by the Master equation in Eq. (§ ex¬ 
cept maybe high-temperature baths. As anticipated in 
Sec. [TTj however, from the perspective of the stability of 
the Majorana mode, one can interpret the dynamics in¬ 
duced by this Master equation with incoherent losses and 
gains as a worst-case scenario. 

We proceed as before by comparing the homogeneous 
case with the disordered case for a non-interacting, i.e. 
[7 = 0, and interacting chain with U ^ 0. We will focus 
on a single protected site at each edge for simplicity. Our 
main results for this setup are summarized in Fig. In 
principle, the overall picture remains identical compared 
to the case with pure losses. Again, we find that the 
presence of disorder suppresses the decoherence process 
of the Majorana mode. At first glance, also the qualita¬ 
tive picture is not changed, as we find a clear exponential 
decay in the clean case and a stretched exponential in the 
disordered case. This is not surprising, as from a sym¬ 
metry point of view, both particle loss alone and particle 
loss and gain break the parity symmetry and in conse¬ 
quence lead to a loss of coherence of the Majorana edge 
mode. 

However, apart from the first observation, that the ad¬ 
ditional particle gain does not completely change the 
overall impact of disorder on the decay of the Majo¬ 
rana, we find subtle modifications in the universality of 
the stretched exponential scaling. In contrast to previ¬ 
ous case with pure losses where the exponent has been 
a 2/3, see Eq. (|^, we find a modified exponent a now 
in the presence of the additional particle gains. Specifi¬ 
cally, from our numerical simulations the exponent turns 
out to be close to a ^ 3/4. However, we have to em¬ 
phasize that the situation is much less clear compared to 
the pure loss case. In particular, we find some slight de¬ 
pendence on the microscopic parameters. Therefore, we 
cannot conclude from our data that the exponent is uni¬ 
versal in this case.l^ For many cases, we find a stretched 
exponential decay with approximately a = 3/4, but this 
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Figure 6. Impact of particle loss and gain on the edge- 
edge correlation 0 in the clean (upper panel: a,b) and dis¬ 
ordered case (lower panel: c,d) Ji G [0,e]: (a) with a fi¬ 
nite external field /x = 0.02,0.03,0.04 {U = 0) and n = 1/2 
for one protected site d = 1; (b) with an on-site repulsion 
U = 0.04, 0.05,0.06 (/i = 0) and = 1 for d = 1. 


value can vary (for the parameters investigated) within 
a range ± 10 %. 

Furthermore, it is worth mentioning that the balanced 
particle loss and gain leads to a more complex picture, 
compared to the pure loss scenario. For example, we 
observed partially decoupled bulk and edge dynamics in 
the case of an interacting chain of length N = 7^ where 
the distance exponential decay in the homogeneous case 
does not depend on the number of protected sites. This 
feature is a result of the uniform bulk dynamics. As 
soon as the bulk settles to a steady state due to gain and 
loss, the absolute length of the bulk is not of importance 
anymore. This stands in contrast to the inhomogeneous 
case, where such a distance dependence is tractable. 

VI. CONCLUSIONS AND OUTLOOK 

We have presented a numerical study of the influence 
of a symmetry-breaking Markovian bath on the stabil¬ 
ity of Majorana edge modes in a Kitaev chain, including 
both disorder and interactions. In a clean setting, the 
occupation of the Majorana mode decays always expo¬ 
nentially as a function of time regardless of the system 
being interacting or not. The corresponding decay rate 
F displays an exponential dependence on the number of 
sites close to the edge where dissipation is not acting, a 
direct consequence of the exponential localization of the 
mode wave-function in the pure loss case. This indicates 
that protecting even a very small number of sites can 
already improve the Majorana mode lifetime by several 
orders of magnitude with respect to the non-protected 
case, if gain is negligible or other interactions prevent an 


































effectively reduced bulk length. 

Strikingly, an even more solid protection mechanism is 
provided by the presence of disorder. In the closed sys¬ 
tem case without dissipation, many-body localization in 
disordered systems pre vents thermalization even in the 
presence of interaction^^^^. Importantly, it has been 
shown that topological properties can be present over 
the ent ire Ha miltonian spectrum in many-body localized 
system^^^^, and that bulk transport is severely inhib¬ 
ited in such phases - which generically display vanishing 
conductivity. We found that the consequences of these 
disorder-induced mechanisms upon the introduction of a 
Markovian bath are that the Majorana mode occupation 
decays in a qualitative different way than in the clean 
case. In the case of loss-only dissipation, which is di¬ 
rectly relevant to cold atom systems, the decay form is a 
stretched exponential, with universal power ^2/3, while 
in the case of both loss and gain - closer to solid state sce¬ 
narios - the decay is still compatible with a stretched ex¬ 
ponential, but with a larger power close to 3/4, whose 
determination is generically more challenging. 

Overall, our results point toward the possibility of us¬ 
ing many-body localization as a stabilizer mechanism of 
edge states in symmetry-protected topological states in 
the presence of a bath with explicitly breaks the symme¬ 
try. While the results presented herein are obtained using 
a controlled numerical approach, it would be of imme¬ 
diate relevance to develop a more general understanding 


using approximate solutions, which capture the main fea¬ 
tures of many-body localized systems coupled to baths, 
e.g., by exploiting the emergence of approximate integrals 
of motion, or by treating the bulk as an effective bath 
for the edges. Moreover, while the work here deals with 
symmetry-protected topological phases, the addition of 
disorder could serve as a stabilizing mechanism under 
dissipation even in phases supporting topolc^ical order 
in combination to long-range entanglemenfP^, where it 
has been recently shown h ow dis order can be used to lo¬ 
calize unwanted excitation d^^ l ^^ l, and could also help in 
further stabiliz ing Maj orana modes against symmetry¬ 
preserving nois d^^ l ^Q ^^. 
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